The decay of magnet ohydro dynamic turbulence in a confined domain 
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The effect of non periodic boundary conditions on decaying two-dimensional magnetohydrody- 
namic turbulence is investigated. We consider a circular domain with no-slip boundary conditions 
for the velocity and where the normal component of the magnetic field vanishes at the wall. Dif- 
ferent flow regimes are obtained by starting from random initial velocity and magnetic fields with 
varying integral quantities. These regimes, equivalent to the ones observed by Ting, Matthaeus and 
Montgomery [Phys. Fluids 29, 3261, (1986)] in periodic domains, are found to subsist in confined 
domains. We examine the effect of solid boundaries on the energy decay and alignment properties. 
The final states are characterized by functional relationships between velocity and magnetic field. 

PACS numbers: 95.30.Qd, 52.65.Kj, 47.ll.Kb 



I. INTRODUCTION 

The influence of initial conditions on decaying mag- 
net ohydro dynamic (MHD) turbulence received consider- 
able interest in the 1980's, because of its relevance to 
explain solar- wind data 0, 0, [I]- Indeed, in magneto- 
hydrodynamics the behavior of decaying turbulent flow 
depends strongly on the initial conditions, and differ- 
ent initial values and ratios of integral quantities can 
lead to a wide variety of distinct behaviors. The first 
systematic study of the different possible types of de- 
cay was performed by Ting, Matthaeus and Montgomery 
[l|, who identified four classes of possible decay be- 
havior, corresponding roughly to a magnetically domi- 
nated, a hydrodynamically dominated, a magnetically- 
hydrodynamically equipartitioned and an erratic transi- 
tion regime. Their study considered the two-dimensional 
case, which is not only relevant in applications in which 
an externally imposed field renders the flows quasi two- 
dimensional, but also from a general physical understand- 
ing of MHD turbulence, which behaves quite similar in 
two and three dimensions, due to the equivalent role of 
the ideal invariants [5]. 

Whereas the influence of the initial conditions on de- 
caying MHD turbulence has been studied and understood 
to some extend, studies on the effect of boundary con- 
ditions have been limited to low resolutions 0, HI, 
imposed by the numerical methods used to account 
for boundaries. Even though these investigations high- 
lighted interesting physics, higher resolution simulations 
are needed to obtain a better understanding of wall- 
bounded MHD, which plays a dominant role in geophys- 
ical flows in the core of planets such as the earth and 
industrial processes involving liquid metals. For the hy- 
drodynamic case it was found that, boundary conditions 
have a significant influence on two-dimensional turbu- 
lence In contrast to the periodic domain, where gen- 
erally a long lasting state is found with a functional sink 
relationship between the vorticity and the stream func- 



tion [T^, corresponding to two counterrotating vortices, 
in a bounded domain with no-slip wall conditions the fi- 
nal state yields an axisymmetric vorticity distribution, 
for which a linear relation between vorticity and stream 
function is observed pU ]. 




27r 

FIG. 1: The computational domain is a square box 2n. The 
fluid domain is a circular container with radius ^ = |§7r, 
surrounded by the solid domain Qg. 

In the present work we propose an extension of the vol- 
ume penalization method [l2[ to two-dimensional MHD 
to compute decaying flows in bounded domains using an 
efficient Fourier pseudo- spectral method. We address the 
following questions: what is the influence of confinement 
by fixed solid boundaries on decaying two-dimensional 
MHD turbulence? Do the four regimes found by Ting 
et al. [l| continue to exist in the presence of boundaries? 
What are the final (viscously decaying) states? 



II. GOVERNING EQUATIONS AND 
NUMERICAL METHOD 



We consider resistive MHD, formulated in usual di- 
mensionless variables u = (u^ v) and B = {Bx^ By) which 
are respectively the velocity and the magnetic field. The 
flow is considered to be two-dimensional, incompressible 
and we assume the mass density to be constant. 
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The governing equations are the fohowing: 
^ + u • Vu = -Vp + j X B + i/V^u - \{u - Uo) (1) 

— = V X (u X B) + T/V^B - -x(B - Bo) (2) 
V-u = V-B = 0(3) 

Here u and r] are respectively the kinematic viscosity 
and the magnetic diffusivity. uoGz = V x u is the vortic- 
ity, jez = V X B is the current density. Furthermore we 
define the vector potential a — aGz ^s B = V X a and 
the stream function ip as u = V^V^ = {—dip/dy^ dip/dx). 
An originality in our approach is the way in which the 
boundary conditions are imposed: we use volume (or sur- 
face in 2D) penalization [12, 13] to include the boundary 
conditions. This method has the advantage that arbi- 
trary basis-functions can be used. In our case a Fourier 
pseudo-spectral code is employed. The advantage with 
respect to a method based on a decompostion in terms 
of Chandrasekhar-Kendall eigenfunctions 0,0,0] is, that 
fast Fourier transforms can be used, allowing for high res- 
olution computations of low computational cost. Also, 
its application to three dimensional flows is conceptually 
straightforward and will be addressed in a future work. 
The additional terms on the right hand side of equa- 
tion ^ and ([2]) correspond to this penalization-method. 
The quantities uo and Bo correspond to the values im- 
posed in the solid part of the numerical domain illus- 
trated in figure [TJ Here we choose uo = and Bo = By 
(where By is the tangential component of B at the wall), 
corresponding to vanishing velocity and no penetration 
of magnetic field into the solid domain which is hence 
considered as a perfect conductor, coated inside with a 
thin layer of insulant, which guarantees that the current 
density cannot penetrate into the solid Q. The mask 
function x is equal to inside Qf (where the penaliza- 
tion terms thereby dissappear) and equal to 1 inside l^g. 
The physical idea is to model the solid part as a porous 
medium whose permeability e tends to zero Bill. For 
e ^ 0, where the obstacle is present, the velocity u tends 
to Uo and the magnetic field B tends to Bo- The nature 
of the boundary condition for the velocity is thus no-slip 
at the wall. 

In the two-dimensional case it is convenient to take the 
curl of ^ and ([2]) to obtain after simplification equa- 
tions for the vorticity and current density. These are 
scalar equations which automatically satisfy the incom- 
pressibility conditions (|3]). The equations are then 

duj 1 

= -u • Vo) + B • Vj + vV'^Lo - -V X [x(u - uo)](4) 

^ = -V^ (u X B) + riV^j --Vx [x(B - Bo)] (5) 

The equations are discretized with a classical Fourier 
pseudo-spectral method imposing periodic boundary 
conditions on the square domain of size 27r, using 512^ 
grid points. At each iteration the fields are dealiased by 



spherical truncation following the 2/3 rule. The penal- 
ization parameter e, corresponding to the permeability of 
the solid domain, is taken equal to 10~^, a value validated 
by a systematic study of the sensitivity of the results 
to this parameter [l3|. The fluid viscosity u and mag- 
netic diffusivity r] were taken equal to 10~^, the timestep 
dt equals to 5.10"^. The initial kinetic and magnet ic 
Reynolds number are defined as Re = 2ry^2Eu{t = 0)/z/ 
and Rcm = 2r^/2EB{t = 0)///, where r is the radius of 
the domain and, Eu and Eb are the kinetic and magnetic 
energies, respectively (see table |I]). 



III. INITIAL CONDITIONS 

Both vorticity and current density fields are initial- 
ized with Gaussian random initial conditions. There- 
fore, their Fourier transforms Q and j, where cD(k) = 
J C(;(x)e~*^"^dx, are initialized with random phases 
and, their amplitudes give the energy spectra: 



Eu(k),EB(k) (X 



k 



{g^{k/ko)y 



(6) 



with /c = |k| and, where g = 0.98 and ko = jV^tt. This 
energy spectrum follows a power law proportional to k~^ 
at large wavenumbers and was chosen to compare with 
simulations performed in the periodic case. Both fields 
are statistically identical. The corresponding fields u and 
B are calculated from co and j using the Biot-Savart law. 

For vanishing viscosity and resistivity, two-dimensional 
MHD has three conserved invariants. The total energy 
is E^ defined as sum of the kinetic energy Eu and the 
magnetic energy Eb ' 



E = Eu -\- Eb = 
He is the cross helicity: 



B|") d'x (7) 



VL-'B d?x 



(8) 



which mesures the global correlations between u and B 
and A is the integral of the squared vector potential: 



A=- i a" d'x 



(9) 



As was shown by Ting et al. [l^ for periodic boundary 
conditions, the dynamics of decaying MHD turbulence 
depend strongly on the initial values of these invariants. 
Because of its interest for the present study we recall 
briefly the four distinct decay regimes discerned by Ting 
et al. [H depending on the initial values and ratios of 
the invariants. First, in the case of small initial He and 
Eb > Eu, di magnetically dominated regime is obtained. 
Selective decay is observed in this regime which corre- 
sponds to the decay of Eu relative to Eb- Second, in 
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tj/A 


Eu/Eb 


He 


Re 


Rem 


Regime I 


UA 


0.068 


0.05 


3580 


7310 


Regime II 


3- 10^ 


2.6 • 10^ 


3.6 • 10"^ 


7900 


53 


Regime III 


53.5 


1.48 


0.23 


5340 


4620 


Regime IV 


14.03 


1.18 


9.5 • 10-^ 


5970 


5340 



TABLE I: Initial values of the four different regimes. 



the case of vanishingly small initial magnetic energy, the 
Lorentz force acting in the vorticity equation can not 
become strong enough, so that the vector potential is 
advected like a passive scalar. Following Biskamp and 
Welter [ij] , the magnetic field may however be amplified 
even if the initial ratio Eb/Eu is very small, given that 
T] is sufficiently small. They found that Eu/Eb < 
is necessary such that the magnetic field can be inten- 
sified. This is a regime which essentially corresponds to 
the Navier-Stokes limit. Third, in the case of substan- 
tial initial cross-helicity, the turbulence tends towards 
an Alfvenic state in which u and B are aligned or anti- 
aligned and approximately equipartitioned. This process 
is called dynamic alignment and the ratio E/\Hc\ tends to 
two. This is a state free from nonlinear interactions, in- 
hibiting cascade processes (even though this depletion of 
nonlinearity is rather slow with increasing cross-helicity 
U). The fourth and final regime is an erratic regime 
which might tend to different final states and which could 
be related to various competing subregions with unequal 
sign of cross-helicity. This regime can be found if the 
flow is initialized with small cross-helicity and compara- 
ble kinetic and magnetic energies. 

Whether these regimes persist in the presence of solid 
boundaries is one of the main questions we want to an- 
swer in the present work. To obtain the desired initial 
conditions corresponding to the four regimes we proceed 
as follows: 

Starting from random initial conditions in Fourier 
space, we renormalize u and B in physical space by vary- 
ing the coefficient a: 
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(10) 



FIG. 2: Time evolution of integral quantities. Top: ratio 
of kinetic and magnetic energy, Eu/Eb- Center: ratio of 
total energy and integral of the squared vector potential, E/A. 
Bottom: ratio of total energy and magnetic helicity, E/\Hc\. 



This generally yields initial conditions with vanishingly 
small cross-helicity, and initial conditions for regime I,II 
and IV can hereby be created. In the case of regime 
III, a non-zero cross-helicity needs to be imposed. We 
achieve this by creating a random initial condition for u 
and a perpendicular field by rotating u by 7r/2. The 
magnetic field is then obtained by a linear combination 
of the 2 fields: 



B* =/3U+(l-;5)U^ 



(11) 



Hereby any given cross helicity can be imposed. Table 
H] summarizes the initial values of E / Eu/Eb and He 
for the four different regimes, together with the Reynolds 
number. 



IV. RESULTS AND DISCUSSION 
A. Characterization of the different decay regimes 

In figure [2] we show the time evolution of several inte- 
gral quantities for the 4 different sets of initial conditions. 
The main observation is that the 4 different regimes, dis- 
cerned by Ting et al. [l| are robust enough to survive 
within a bounded domain. We now discuss the results in 
more detail. 

Regarding the ratio of kinetic and magnetic energy 
(Fig. [21 top), it is observed that in the absence of initial 
cross-helicity (case I, II and IV) the magnetic energy fi- 
nally dominates, unless it is very small initially (Navier 
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FIG. 3: Time evolution of Eu/ Eb in a periodic domain, start- 
ing from similar initial conditions as in figure [2l 



Stokes limit). However, if the initial cross-helicity is ini- 
tially large and Eu/Eb is of order unity, the flow energy 
will remain approximately equipartitioned between the 
velocity and magnetic field. 

This picture is confirmed by the time evolution of the 
ratio E/A (Fig. O center). In this representation it is 
however emphasized that in the Navier- Stokes limit (case 
II), the character of the magnetic field has changed: in 
the ideal system (vanishing viscosity and magnetic diffu- 
sivity), A is a quantity that cascades towards the small 
wavenumbers. In a non-ideal system an inverse cascade 
generally slows down the dissipation rate of the quan- 
tity. However, in the limit of small Lorentz force, the 
equations of the vorticity and vector potential become 
equivalent to the equations that describe a passive scalar 
advected by a two-dimensional velocity field. The passive 
scalar being a quantity which cascades towards higher 
wavenumbers, the vector potential gets dissipated faster 
in this case than in the case where the Lorentz force is 
significant. This results in a rapid increase of the quan- 
tity E/A in case II. 

The ratio E/\Hc\ (Fig. O bottom) attains its minimum 
absolute value 2 for case HI. This corresponds to dynamic 
alignment: the velocity field is equal in magnitude and 
perfectly aligned, or anti-aligned with the magnetic field. 
The erratic regime is clearly represented by case IV, in 
which the cross-helicity approaches a value close to zero. 
As we will see in the following, this is caused by different 
subregions with oppositely valued He. 

For comparison, we show in figure [3] in a pe- 

riodic domain, starting from similar initial conditions as 
in the bounded case using the same numerical param- 
eters. Even though the trend is similar, we see that a 
more oscillatory behavior for case I and II is observed 
than in the case of the bounded domain. This oscilla- 
tory behavior is related to energy exchange between the 
magnetic field and the velocity field by means of Alfven 
waves Whereas in a periodic domain these waves 

can freely propagate, in a bounded domain they might 
be more rapidly suppressed, explaining the less oscilla- 
tory behavior of Eu/Eb in a bounded domain. Further 



FIG. 4: Time evolution of the average alignment cos 0, be- 
tween the magnetic field and the velocity field. 
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FIG. 5: Probability density of cos 6* at t = 40, t = 450 and 
t — 1250 in the regime IV. 



research is needed to clarify this. 

The quantity E/Hc gives a measure for the dynamic 
alignment, which corresponds to measuring both the 
equipartitioning of energy and the alignment properties. 
If we are exclusively interested in the alignment proper- 
ties, the relative cross helicity, which corresponds to the 
cosine of the angle between the velocity and magnetic 
field vector, 



cos 



[EuEbY^^ 



(12) 



should be considered. 

In figure m cos6> is plotted as a function of time. It can 
be observed that in case I and HI, the velocity field tends 
to a nearly aligned state. In case II and IV, this quantity 
remains close to zero, however for a different reason. In 
case II, the alignment is small, because the vector po- 
tential is advected as a nearly passive scalar. In case IV 
the local alignment is large but different aligned or anti- 
aligned regions cancel out the contributions, yielding a 
net-global alignment close to 0. This can be observed 
in the corresponding probability distribution function of 
cosO at t = 40 and t = 450, shown in figure [5l Nev- 
ertheless, for long time {t = 1250) we observe an anti- 
alignment. 
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FIG. 6: Time evolution of the total energy in log- log scale 
(top) and in log-lin scale (bottom). The solid line (top) cor- 
responds to t~°'^ and the dotted line (top) corresponds to 



B. Energy decay and visualizations 

The decay of total energy is shown in figure [6l At 
intermediate times, the energy in case I and II decays 
following a powerlaw with exponents varying for the dif- 
ferent sets of initial conditions (Fig. [6l top). The expo- 
nents of these powerlaws are approximately —0.6 (dotted 
line) for regime I and —0.4 (solid line) for regime IV. It 
is seen that these powerlaws are observed only after an 
initial period of rapid decay. In the other cases no clear 
powerlaw behavior can be identified. This can be com- 
pared to previous studies [sl, [IB] in which values around 
—0.75 and —1 were found for the decaying periodic case. 
In case III, in which dynamic alignment is observed, no 
clear power-law behavior is observed. In this case the 
nonlinear interactions are progressively damped by the 
alignment process, so that no selfsimilar period is ob- 
served in the energy decay. At late times (Fig. [6l bot- 
tom) all cases show an exponential viscous decay of the 
form E ~ g-2azyt with a = 1.5 in case I and a = 2 in 
cases II, III and IV, a value related to the largest Stokes 
eigenmode of the circle {a = 1.64), which contains most 
of the energy, as found in [ll| for the hydrodynamical 
case. 

Figures [3 and [8] show the vorticity and the current den- 
sity field, respectively. For each of the cases I-IV, three 
typical time instants are visualized. These instants are 



t = 5, showing the self-organization of the flow at early 
times, t = 40, when nonlinear processes are dominating 
and t = 250 (regime I) and t = 1250 (regimes II, III and 
IV), corresponding to the final, viscously decaying state. 

One flagrant feature of the visualizations is the local 
alignment of the magnetic and velocity field. Indeed in 
most regimes the vorticity and current density fields are 
rather similar. We also observe the coincidence of the 
maxima of co and of j which may have some effect on the 
stabilization of vorticity and current filaments. In case 
I an almost perfect axi-symmetrical state is achieved at 
t = 250. Case II is the only case in which the forma- 
tion of circular vortices is well pronounced, leading to a 
roll up of the current sheets. Apparently in the other 
regimes the Lorentz force suppresses the generation of 
circular vortices. Case III shows almost identical mag- 
netic and velocity fields, as expected in this case of dy- 
namic alignment, in which u and B are aligned (or anti- 
aligned) and in which kinetic and magnetic energies are in 
equipartition. Case IV is a typical example of the erratic 
regime: at the intermediate time, four dominant flow 
stuctures are observed, with both positive and negative 
cross-helicity. Locally the flow is close to an aligned or 
anti-aligned state, but globally the cross-helicity is weak 
because the different regions with opposite contributions 
cancel each other out. 



C. Final states 

A supplementary information on the final states is 
given by scatter-plots. It was shown by Joyce and 
Montgomery [l^ that in hydrodynamic unbounded two- 
dimensional flows a long lasting final state is reached, 
depleted from nonlinearity. This state is characterized 
by a functional relation between the vorticity and the 
streamfunction of the form co ~ sinh(?/^). That a func- 
tional relation leads to a state, depleted from nonlinearity 
is easily shown from the equation for the vorticity: 

{dt-J^A)uj = [u,i^], (13) 

with the Poisson bracket defined as [a, 6] = 
{da/dx){db/dy) — {da/dy){db/dx). A functional 
relation co = F(?/;) leads to a vanishing Poisson bracket. 
If we consider now the equations for incompressible 
MHD: 

{dt-iyA)co = [u,^]-[aj] (14) 
{dt-vA)a = [a,V^], (15) 

we see that two nonlinear it ies play a role: [a;, ip] and 
[a,j]. The term [a^tjj] can be considered as a pseudo- 
nonlinearity if tjj is regarded as given. Although impor- 
tant theoretical progress has been made in the compre- 
hension of final states [13] no analytical nontrivial solu- 
tion is presently known for the case of decaying MHD 
turbulence. It was however shown in Kinney et al. [l^ 
that close to functional relations do exist in homogeneous 



FIG. 7: Vorticity at different instants in the circular domain. From top to bottom: regime I, regime II, regime III and regime IV; from 
left to right: t = 5, t = 40 and in the last column the time coresponds to t = 250 for regime I and t = 1250 for regimes II, III and IV. 



two-dimensional MHD turbulence. In figure [9] we show 
for the cases I-IV these scatter plots corresponding to the 
three nonlinearities. 

In case I we see a well defined nonlinear functional rela- 
tion u;{tlj). Clearly, we have a non trivial final state. The 
plot a vs. j shows a straight line, which corresponds to a 
vanishing Lorentz-force: the magnetic field does not in- 
teract with the velocity field at this final period of decay. 
The plot a vs. ip also shows a clear functional relation. 
In case II, the scatter plots do not show such clear func- 
tional relations which is due to the fact that the fiow is 
not yet sufficiently relaxed. The plot uj vs. i/j is per- 
haps closest to a functional relation. In case III we see 
as expected a vanishing nonlinearity: for dynamic align- 
ment it can be shown that nonlinearities vanish in the 



perfectly aligned case, when the equations are stated in 
Elsasser variables (see for example [4]). In case IV it is 
expected that eventually the same behavior is observed 
as in case I. If the initial Reynolds number is initially too 
low this behavior will however not be observed. Prelim- 
inary computations were performed at lower resolution, 
which showed that non-trivial final states are only ob- 
served if the initial Reynolds number is sufficiently high. 
Otherwise linear relations are obtained for all different 
scatter plots. 

V. CONCLUSION 

We have investigated the inffuence of non-periodic 
boundary conditions on decaying two-dimensional mag- 



FIG. 8: Current density at different instants in the circular domain. From top to bottom: regime I, regime II, regime III and regime IV; 
from left to right: t = 5, ^ = 40 and in the last column the time coresponds to t = 250 for regime I and t = 1250 for regimes II, III and IV. 



net ohydro dynamic turbulence. The use of a penalization 
method in combination with a classical Fourier pseudo- 
spectral method allows for efficient resolution of MHD 
flows in bounded domains. 

A main result is the observation of the robustness of 
the four different regimes discerned by Ting et al 
The same trends are found as in their pioneering work, 
depending on the initial values of the kinetic energy, 
magnetic energy, vector potential and cross-helicity. A 
detailed description was given of the relaxation-process 
which leads to the final states. In the case of a magneti- 
cally dominant, cross-helicity free case, a clear nontrivial 



functional relation was observed describing the magnetic 
and velocity fields. Functional relationships were also 
observed in regimes III and IV, while in regime II this 
functional relation was less clear. 

Future work will address the influence of other types 
of boundary conditions for the magnetic field and also 
other geometries will be studied. 
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